Regresión Lineal Múltiple
TP Final
Alumnos
Nicolás Dominutti, Carlos Suarez Gurruchaga, Hernán Telechea
Primer acercamiento y correcciones
SegVial <- read_excel("SegVial.xlsx")
(head(SegVial,3))## # A tibble: 3 × 20
## `Codigo Pais` Pais Ciudad Poblacion DenPob ArCiclista ArBajaVel PMPeatones
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 ES Spain Barcelo… 1620343 16014 0.15 0.775 32.0
## 2 ES Spain Madrid 3223334 5332 0.067 0.065 34.0
## 3 FR France Bordeaux 254436 5104 0.106 0.239 21.0
## # … with 12 more variables: PMCiclistas <chr>, PMTPublico <chr>,
## # PMVMotor <chr>, PeatAuto <dbl>, CicAuto <dbl>, V2RMSM <dbl>,
## # V2RMAuto <dbl>, AutoSM <dbl>, AutoAuto <dbl>, Precipitacion <chr>,
## # Temp <chr>, PBI <dbl>glimpse(SegVial)## Rows: 24
## Columns: 20
## $ `Codigo Pais` <chr> "ES", "ES", "FR", "FR", "FR", "FR", "FR", "FR", "FR", "F…
## $ Pais <chr> "Spain", "Spain", "France", "France", "France", "France"…
## $ Ciudad <chr> "Barcelona", "Madrid", "Bordeaux", "Lille", "Lyon", "Mar…
## $ Poblacion <dbl> 1620343, 3223334, 254436, 232787, 516092, 863310, 285121…
## $ DenPob <dbl> 16014, 5332, 5104, 6687, 10758, 3564, 4997, 4704, 4586, …
## $ ArCiclista <chr> "0.15", "0.067", "0.106", "0.098", "0.109", "0.044", "0.…
## $ ArBajaVel <chr> "0.775", "0.065", "0.239", "0.634", "0.284", "0.118", "0…
## $ PMPeatones <chr> "32.0", "34.0", "21.0", "32.0", "34.0", "34.0", "26.0", …
## $ PMCiclistas <chr> "2.0", "0.46", "3.0", "2.0", "2.0", "1.0", "2.0", "5.0",…
## $ PMTPublico <chr> "39.0", "24.35", "9.0", "10.0", "19.0", "11.0", "8.0", "…
## $ PMVMotor <chr> "27.0", "39.89", "67.0", "56.0", "45.0", "54.0", "64.0",…
## $ PeatAuto <dbl> 22, 201, 13, 5, 31, 56, 21, 13, 20, 58, 5, 12, 153, 60, …
## $ CicAuto <dbl> 5, 18, 5, 1, 11, 9, 4, 7, 2, 14, 8, 2, 54, 14, 28, 15, 2…
## $ V2RMSM <dbl> 29, 95, 2, 2, 12, 45, 8, 6, 6, 48, 4, 5, 6, 8, 1, 1, 2, …
## $ V2RMAuto <dbl> 96, 275, 14, 7, 39, 146, 21, 13, 50, 132, 8, 29, 55, 11,…
## $ AutoSM <dbl> 0, 44, 3, 2, 8, 30, 9, 1, 5, 12, 4, 5, 20, 18, 6, 7, 11,…
## $ AutoAuto <dbl> 4, 72, 1, 6, 21, 48, 12, 6, 6, 24, 3, 12, 122, 31, 11, 1…
## $ Precipitacion <chr> "565.0", "420.9", "944.1", "742.5", "831.9", "515.4", "6…
## $ Temp <chr> "18.2", "15.0", "13.8", "10.8", "12.5", "15.5", "15.1", …
## $ PBI <dbl> 31070.56, 35290.16, 35464.54, 31291.93, 48140.02, 36231.…Tenemos 20 variables analizadas con 24 datos cada una, no se observan valores NaNs. 11 de esas variables aparecen como categóricas, el resto como cuantitativas. Sin embargo, al observar la tabla de datos, solo 3 deberían ser categóricas. Las demás aparecen como tal por tener decimales. Las vamos a transformar.
SegVial[,c(6, 7, 8, 9, 10, 11, 18, 19)] = lapply(SegVial[,c(6, 7, 8, 9, 10, 11, 18, 19)], function(x) as.numeric(x))
summary(SegVial)## Codigo Pais Pais Ciudad Poblacion
## Length:24 Length:24 Length:24 Min. : 232787
## Class :character Class :character Class :character 1st Qu.: 432558
## Mode :character Mode :character Mode :character Median : 542400
## Mean : 976996
## 3rd Qu.: 932826
## Max. :3600203
## DenPob ArCiclista ArBajaVel PMPeatones
## Min. : 1417 Min. :0.02700 Min. :0.0240 Min. : 4.00
## 1st Qu.: 3223 1st Qu.:0.06425 1st Qu.:0.1003 1st Qu.:11.50
## Median : 4212 Median :0.10200 Median :0.2380 Median :26.50
## Mean : 5483 Mean :0.11675 Mean :0.2974 Mean :23.98
## 3rd Qu.: 5161 3rd Qu.:0.15325 3rd Qu.:0.4098 3rd Qu.:33.25
## Max. :20772 Max. :0.33300 Max. :0.8710 Max. :47.00
## PMCiclistas PMTPublico PMVMotor PeatAuto
## Min. : 0.460 Min. : 8.00 Min. :17.00 Min. : 5.00
## 1st Qu.: 1.000 1st Qu.:11.75 1st Qu.:43.72 1st Qu.: 18.25
## Median : 2.000 Median :19.00 Median :53.00 Median : 36.50
## Mean : 2.811 Mean :21.84 Mean :51.39 Mean : 70.71
## 3rd Qu.: 3.000 3rd Qu.:31.50 3rd Qu.:63.25 3rd Qu.: 67.25
## Max. :14.000 Max. :39.00 Max. :71.00 Max. :398.00
## CicAuto V2RMSM V2RMAuto AutoSM
## Min. : 1.00 Min. : 1.00 Min. : 2.00 Min. : 0.00
## 1st Qu.: 5.75 1st Qu.: 3.50 1st Qu.: 15.50 1st Qu.: 4.75
## Median : 14.00 Median : 6.00 Median : 29.50 Median : 7.50
## Mean : 30.46 Mean : 21.96 Mean : 74.83 Mean : 18.04
## 3rd Qu.: 27.25 3rd Qu.: 20.00 3rd Qu.: 65.25 3rd Qu.: 21.00
## Max. :340.00 Max. :143.00 Max. :371.00 Max. :138.00
## AutoAuto Precipitacion Temp PBI
## Min. : 1.00 Min. : 420.9 Min. : 6.800 Min. : 23867
## 1st Qu.: 6.00 1st Qu.: 637.9 1st Qu.: 9.925 1st Qu.: 34399
## Median : 17.50 Median : 800.3 Median :11.100 Median : 37087
## Mean : 40.17 Mean : 753.5 Mean :12.021 Mean : 50517
## 3rd Qu.: 49.25 3rd Qu.: 829.6 3rd Qu.:14.475 3rd Qu.: 48993
## Max. :217.00 Max. :1245.3 Max. :18.200 Max. :152178Estudio de las covariables relacionadas con la poblacion
Vamosa empezar nuestro estudio de las covariables poniendo foco en la población de cada ciudad
boxplot(SegVial$Poblacion, main = "Poblacion", col = "blue", horizontal = T)
A priori, vemos cómo existen algunas ciudades que sobresalen en cuanto a población se refiere.
SegVial %>% ggplot(aes(x=reorder(Ciudad,Poblacion,max), y=Poblacion))+ geom_bar(stat = 'identity', col='red') +
geom_bar(stat = 'identity', aes(x=Ciudad, y=PBI), col='blue') +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="Población y PBI en cada ciudad")+
xlab('Ciudades')
En este gráfico podemos ver la Población de cada ciudad en rojo y el PBI de las mismas en azul. Observamos que las 4 ciudades de las que hablábamos son: Inner London, Madrid, Rome y Paris, teniendo Inner London y París 2 de los PBI más elevados. Quizá el caso que más llama la atención es Montpellier, que presenta baja población pero muy alto PBI.
aux_plot = copy(SegVial)
aux_plot$PBI_per_capita = aux_plot$PBI / aux_plot$Poblacion
aux_plot %>% ggplot(aes(x=reorder(Ciudad,PBI_per_capita,max), y=PBI_per_capita))+ geom_bar(stat = 'identity', col='green') +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="PBI per capita en cada ciudad")+
xlab('Ciudades')
Vemos como al calcular el PBI per cápita Montpellier destaca y Madrid pasa a ser la ciudad con el indicador más bajo.
Analicemos ahora la Densidad Poblacional de cada ciudad
boxplot(SegVial$DenPob, main = "Densidad Poblacional", col = "blue", horizontal = T)
SegVial %>% ggplot(aes(x=reorder(Ciudad,DenPob,max), y=DenPob))+ geom_bar(stat = 'identity', col='red') +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="Densidad Poblacional de cada ciudad")+
xlab('Ciudades')
Observamos 4 ciudades que presentan valores de densidad poblacional mas elevadas que el resto: Paris, Barcelona, Inner London y Lyon. Sabemos que ciudades como Inner London o París tienen una elevada población, pero ciudades como Barcelona no, pero por otro lado también sabemos que la densidad poblacional es un ratio entre población sobre km2, veamos cómo son los km2 para cada ciudad.
aux_plot$km2 = aux_plot$Poblacion / aux_plot$DenPob
aux_plot %>% ggplot(aes(x=reorder(Ciudad,km2,max), y=km2))+ geom_bar(stat = 'identity', col='black') +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="km2 de cada ciudad")+
xlab('Ciudades')
Como decíamos antes, las ciudades con menor área van a ser más propensas a tener una concentración poblacional mayor, en este caso vemos cómo Barcelona, Lyon y Lillie tienen alta densidad poblacional no por su cuantiosa población, sino por su reducida área de extensión.
Estudio de las variables target
Este trabajo tiene como objetivo intentar darle explicación a las variables de los accidentes ocurridos, con el fin de proponer políticas públicas para disminuirlos a futuro. En este apartado vamos a investigar un poco estas features target.
Primero que nada vamos a analizar cómo son las incidencias de cada tipo de accidente en todo el dataset y luego por ciudad
color <- brewer.pal(6, "Pastel1")
pie(c(sum(SegVial$PeatAuto), sum(SegVial$CicAuto), sum(SegVial$V2RMSM),
sum(SegVial$V2RMAuto), sum(SegVial$AutoSM), sum(SegVial$AutoAuto)),
labels = c(" PEATON->AUTO (27.6%)"," CICLISTA->AUTO (11.9%)", " DOS RUEDAS(M)->SI MISMOS (8.6%)",
" DOS RUEDAS(M)->AUTO (29.2%)", " AUTO->SI MISMOS (7%)", " AUTO->AUTO (15.7%)"),
col = color, main = "DISTRIBUCION DE LOS ACCIDENTES EN EL DATASET")
Tenemos un dataset en donde la mayoría de los accidentes son de: * vehículos de 2 ruedas vs autos * autos a peatones
incidencias = copy(SegVial)
incidencias = incidencias %>% dplyr::select(Ciudad, V2RMSM, V2RMAuto, AutoAuto, AutoSM, PeatAuto, CicAuto, accidentes_viales)
incidencias$V2RMSM = round(incidencias$V2RMSM / incidencias$accidentes_viales,4)
incidencias$V2RMAuto = round(incidencias$V2RMAuto / incidencias$accidentes_viales,4)
incidencias$AutoAuto = round(incidencias$AutoAuto / incidencias$accidentes_viales,4)
incidencias$AutoSM = round(incidencias$AutoSM / incidencias$accidentes_viales,4)
incidencias$PeatAuto = round(incidencias$PeatAuto / incidencias$accidentes_viales,4)
incidencias$CicAuto = round(incidencias$CicAuto / incidencias$accidentes_viales,4)
incidencias = melt(as.data.table(incidencias))## Warning in melt.data.table(as.data.table(incidencias)): id.vars and measure.vars
## are internally guessed when both are 'NULL'. All non-numeric/integer/
## logical type columns are considered id.vars, which in this case are columns
## [Ciudad, ...]. Consider providing at least one of 'id' or 'measure' vars in
## future.incidencias[incidencias$variable!='accidentes_viales',] %>% ggplot(aes(x=Ciudad ,y=value, fill = variable))+ geom_bar(stat = "identity") +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="Incidencia del tipo de accidentes por ciudad")
De este grafico podemos sacar información interesante:
- En Barcelona más del 75% de los accidentes involucran un vehículo de 2 ruedas motorizado
- Madrid, Niza, Marsella, Paris y Toulouse, son los otros países en donde más del 50% de los accidentes involucran este tipo de vehículos
- Los accidentes que involucran a ciclistas representan un porcentaje mayor del total en Bristol, Inner London, Manchester, Oslo y Strasbourg, siendo en Bristol los de mayor incidencia
- Los choques entre autos son preponderantes en Sheffield
- Los accidentes que involucran peatones llegan al 48% en Glasgow, seguido por el 46% en Liverpool y el 42% en Bradford
Ahora vamos a ver más en detalle los valores que toman las variables respuesta en el dataset
par(mfrow = c(3, 2))
boxplot(SegVial$PeatAuto, col = "blue",
main = "PEATON->AUTO", horizontal = T)
boxplot(SegVial$CicAuto, col = "blue",
main = "CICLISTA->AUTO", horizontal = T)
boxplot(SegVial$V2RMSM, col = "blue",
main = "MOTO->SI MISMO", horizontal = T)
boxplot(SegVial$V2RMAuto, col = "blue",
main = "MOTO->AUTO", horizontal = T)
boxplot(SegVial$AutoSM, col = "blue",
main = "AUTO->SI MISMO", horizontal = T)
boxplot(SegVial$AutoAuto, col = "blue",
main = "AUTO->AUTO", horizontal = T)
aux_para_plot = reshape2::melt(SegVial[,c('Ciudad','V2RMSM','V2RMAuto','AutoAuto','AutoSM','PeatAuto','CicAuto','accidentes_viales')], value.name = "accidentes")## Using Ciudad as id variablesaux_para_plot %>% ggplot(aes(x=Ciudad ,y=accidentes, alpha=.2, color = variable))+ geom_point() +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="Cantidad de accidentes por tipo en cada ciudad")
De estos gráfico podemos ver que: * Hay ciudades con cantidades atípicamente altas de accidentes de los diversos tipos * Inner London resulta ser la ciudad que tiene mayor cantidad de accidentes totales, accidentes PeatAuto y accidentes V2RMAuto * Roma y Madrid ocupan el 2do y 3er puesto respectivamente en cuando a accidentes totales
Recordemos del apartado de población, que estas son las 3 ciudades que mayor población tienen.
Veamos ahora con un poco más de detalle por tipo de accidente
Como dijimos, observamos que tanto la ciudad de Inner London, como Roma y Madrid aparecen en la mayoria de los casos con valores muy altos de accidentes de distintos tipos.
Estudio de variables regresoras
Areas ciclista y baja velocidad
Vamos a continuar analizando las covariables relacionadas al area destinada para cada uso
par(mfrow = c(2,1))
boxplot(SegVial$ArCiclista, main = "Area Ciclista vs Motor", col = "blue",horizontal = T)
boxplot(SegVial$ArBajaVel, main = "Area baja velocidad", col = "blue",horizontal = T)
Sobre la relacion entre el area destinada de baja velocidad, no observamos valores alejados del grupo para alguna ciudad en particular.
round(mean(SegVial$ArBajaVel), 2)## [1] 0.3Observamos que en promedio 7 de cada 10 km2 no tienen restricciones de velocidad en las Ciudades.
SegVial %>%
dplyr::select(Pais, Ciudad, ArCiclista) %>%
arrange(-ArCiclista) %>% head(2)## # A tibble: 2 × 3
## Pais Ciudad ArCiclista
## <chr> <chr> <dbl>
## 1 France Strasbourg 0.333
## 2 France Nantes 0.327Respecto a la relacion entre el area destinada para Ciclistas vs Vehiculos Motorizados, observamos 2 Ciudades de Francia, que presentan una cantidad de superficie mayor para el uso de la bicicleta, respecto al grupo de ciudades estudiado en general.
Parcitipación modal
aux_para_plot = reshape2::melt(SegVial[,c('Ciudad','PMPeatones','PMCiclistas','PMTPublico','PMVMotor')], value.name = "PM")## Using Ciudad as id variablesaux_para_plot %>% ggplot(aes(x=Ciudad ,y=PM, alpha=.2, fill = variable))+ geom_bar(stat = "identity") +
theme(axis.text.x=element_text(angle = 90)) +
labs(title="PM por tipo de movilidad en cada ciudad")
A simple vista podemos ver que: * Existe una clara preponderancia del transporte via vehiculo motorizado * París es la ciudad con una menor PM de vehículos a motor y una mayor incidencia de personas que se desplazan caminando. Recordemos que París es una de las ciudades con mayor cantidad de accidentes de V2RM * Barcelona también tiene poco PM de vehículos con motor, preponderando el transporte público * Leeds, Bordeaux, Roma y Sheffield tienen muy bajo PM de peatones * Bristol y Strasbourg tienen niveles elevados de ciclistas, respecto a las demas ciudades
Variables meteorológicas
par(mfrow = c(2,1))
boxplot(SegVial$Temp, main = "TEMPERATURA MEDIA [C]", col = "blue", horizontal = T)
boxplot(SegVial$Precipitacion, main = "PRECIPITACION MEDIA [mm]", col = "blue", horizontal = T)
SegVial %>%
dplyr::select(Pais, Ciudad, Precipitacion) %>%
arrange(-Precipitacion) %>% head(1)## # A tibble: 1 × 3
## Pais Ciudad Precipitacion
## <chr> <chr> <dbl>
## 1 United Kingdom Glasgow 1245.Observamos un valor elevado de mm de agua para la ciudad de Glasgow, ¿Esto impacta en la cantidad de ciclistas?
par(mfrow = c(2, 2))
plot(y = SegVial$PMCiclistas, x = SegVial$Precipitacion, col = "blue",
main = "modal CICLISTAS vs mm PRECIPITADOS",
xlab = "mm PRECIPITADOS", ylab = "CICLISTAS [%]")
plot(y = SegVial$PMPeatones, x = SegVial$Precipitacion, col = "blue",
main = "modal PEATONES vs mm PRECIPITADOS",
xlab = "mm PRECIPITADOS", ylab = "PEATONES [%]")
plot(y = SegVial$PMCiclistas, x = SegVial$Temp, col = "blue",
main = "modal CICLISTAS vs Temp Media")
plot(y = SegVial$PMPeatones, x = SegVial$Temp, col = "blue",
main = "modal PEATONES vs Temp Media")
cor(SegVial$PMCiclistas, SegVial[, c(18:19)], method = "spearman")## Precipitacion Temp
## [1,] -0.04857822 -0.08844308No se observa una correlacion entre una ciudad con muchas lluvias o con mucha temperatura y el optar o no por el transporte en bicicleta.
Revisemos lo mismo con peatones
cor(SegVial$PMPeatones, SegVial[, c(18:19)], method = "spearman")## Precipitacion Temp
## [1,] -0.5121323 0.345434Aqui si se observa una correlacion un poco más fuerte. Para las ciudades que tienen una media de temperatura mas elevada, la gente utiliza mas el transporte a pie, lo cual tiene sentido. A su vez, tambien observamos que ciudades que presentan mas cantidad de lluvias, ven reducido el uso del transporte a pie.
cor(SegVial$PMTPublico, SegVial[, c(18:19)], method = "spearman")## Precipitacion Temp
## [1,] 0.1179801 -0.3302114cor(SegVial$PMVMotor, SegVial[, c(18:19)], method = "spearman")## Precipitacion Temp
## [1,] 0.3628454 -0.03156291Observamos que ciudades con mayor cantidad de precipitaciones, y temperaturas medias mas bajas, tienen un aumento del transporte motorizado.
Continuamos el analisis respondiendo algunas preguntas típicas de los accidentes viales:
¿Mayor densidad de poblacion favorece la cantidad de accidentes en general?
par(mfrow = c(2,3))
plot(SegVial$DenPob, SegVial$AutoAuto, main='DenPob vs accidentes Auto-Auto')
plot(SegVial$DenPob, SegVial$AutoSM, main='DenPob vs accidentes AutoSM')
plot(SegVial$DenPob, SegVial$V2RMAuto, main='DenPob vs accidentes V2RM-Auto')
plot(SegVial$DenPob, SegVial$V2RMSM, main='DenPob vs accidentes V2RMSM')
plot(SegVial$DenPob, SegVial$PeatAuto, main='DenPob vs accidentes Peat-Auto')
plot(SegVial$DenPob, SegVial$CicAuto, main='DenPob vs accidentes Cic-Auto')
Por sentido común, esperaríamos que poblaciones con mayor densidad sean propensas a más accidentes, de cualquier tipo, ya que la concentración de gente en zonas más pequeñas deberia impactar en el ordenamiento del tránsito. Sin embargo, los plots no indican esto, por lo que DenPob no pareciera ser una variable muy indicadora del nivel de accidentes.
¿Mayor cantidad de poblacion favorece a la cantidad de accidentes?
par(mfrow = c(2,3))
plot(SegVial$Poblacion, SegVial$AutoAuto, main='Poblacion vs accidentes Auto-Auto')
plot(SegVial$Poblacion, SegVial$AutoSM, main='Poblacion vs accidentes AutoSM')
plot(SegVial$Poblacion, SegVial$V2RMAuto, main='Poblacion vs accidentes V2RM-Auto')
plot(SegVial$Poblacion, SegVial$V2RMSM, main='Poblacion vs accidentes V2RMSM')
plot(SegVial$Poblacion, SegVial$PeatAuto, main='Poblacion vs accidentes Peat-Auto')
plot(SegVial$Poblacion, SegVial$CicAuto, main='Poblacion vs accidentes Cic-Auto')
t(cor(SegVial$Poblacion, SegVial[, c(12:17)],
method = "spearman"))## [,1]
## PeatAuto 0.7598957
## CicAuto 0.5587469
## V2RMSM 0.6363902
## V2RMAuto 0.7401609
## AutoSM 0.6470724
## AutoAuto 0.6930785Observamos que la población tiene relación aproximadamente lineal respecto a los distintos tipos de accidentes. No sabemos si las variables de accidente son realmente índices o totales de accidentes. Si fueran lo segundo, tiene sentido el efecto - donde hay más gente, podés reportar más casos-. Sin embargo, cabe destacar como vimos en los gráficos previos, que esa gente esté más o menos distante, impacta poco.
nota: Cabe aclarar que se evidencia visualmente la presencia de valores atipicamente altos, que fueron analizados previamente en la seccion “ANALISIS DE LOS TARGETS”.
Vamos a analizar si alguna covariable presenta una correlacion lineal contra la variable target_combinacion (sumatoria de todos los accidentes por ciudad). Si fuera asi, esto indicaria que tiene algo que explicar con respecto a alguno de los targets.
df_pairs_1 = SegVial[,-c(1:3, 5, 12:20)]
ggpairs(df_pairs_1, progress = F)
df_pairs_2 = SegVial[,c(18:21)]
ggpairs(df_pairs_2, progress = F)
Observamos una gran correlacion entre nuestro target_combination y poblacion, asi como tambien con PMTPublico y a continuacion lo siguen ArCiclista, Precipitacion y PMCiclistas. Estas covariables probablemente puedan explicarnos bastante de nuestros targets al momento de armar los modelos.
SegVial$PeatAuto_std = (SegVial$PeatAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$CicAuto_std = (SegVial$CicAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$V2RMSM_std = (SegVial$V2RMSM / SegVial$Poblacion) * (10 ^ 6)
SegVial$V2RMAuto_std = (SegVial$V2RMAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$AutoSM_std = (SegVial$AutoSM / SegVial$Poblacion) * (10 ^ 6)
SegVial$AutoAuto_std = (SegVial$AutoAuto / SegVial$Poblacion) * (10 ^ 6)
SegVial$accidentes_viales_std = (SegVial$accidentes_viales / SegVial$Poblacion) * (10 ^ 6)
SegVial[,-c(1:3, 12:17, 21:26)] = as.data.frame(scale(SegVial[,-c(1:3, 12:17, 21:26)]))
head(SegVial)## # A tibble: 6 × 28
## `Codigo Pais` Pais Ciudad Poblacion DenPob ArCiclista ArBajaVel PMPeatones
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ES Spain Barcel… 0.654 2.23 0.404 1.90 0.662
## 2 ES Spain Madrid 2.28 -0.0320 -0.605 -0.924 0.827
## 3 FR France Bordea… -0.734 -0.0802 -0.131 -0.232 -0.246
## 4 FR France Lille -0.756 0.255 -0.228 1.34 0.662
## 5 FR France Lyon -0.468 1.12 -0.0942 -0.0533 0.827
## 6 FR France Marsei… -0.116 -0.406 -0.884 -0.713 0.827
## # … with 20 more variables: PMCiclistas <dbl>, PMTPublico <dbl>,
## # PMVMotor <dbl>, PeatAuto <dbl>, CicAuto <dbl>, V2RMSM <dbl>,
## # V2RMAuto <dbl>, AutoSM <dbl>, AutoAuto <dbl>, Precipitacion <dbl>,
## # Temp <dbl>, PBI <dbl>, accidentes_viales <dbl>, PeatAuto_std <dbl>,
## # CicAuto_std <dbl>, V2RMSM_std <dbl>, V2RMAuto_std <dbl>, AutoSM_std <dbl>,
## # AutoAuto_std <dbl>, accidentes_viales_std <dbl>Como ya dijimos, el objetivo del presente trabajo es encontrar caracteristicas urbanas para encontrar las vulnerabilidades de los usuarios viales, y de esta manera proponer politicas publicas eficaces para disminuir la fatalidad y cantidad de accidentes. Dicho esto, dado las características del dataset, se evidencia que las covariables Pais, Codigo Pais y Ciudad, no nos darian informacion util para este analisis, y por tanto optamos por descartarlas.
df_peatauto = SegVial[, -c(1:3, 12:17, 21, 23:28)]
df_cicauto = SegVial[, -c(1:3, 12:17, 21, 22, 24:28)]
df_v2rmsm = SegVial[, -c(1:3, 12:17, 21, 23, 25:28)]
df_v2rmauto = SegVial[, -c(1:3 ,12:17, 21:24, 26:28)]
df_autosm = SegVial[, -c(1:3 ,12:17, 21:25, 27:28)]
df_autoauto = SegVial[, -c(1:3 ,12:17, 21:26, 28)]
df_accidentesviales = SegVial[, -c(1:3 ,12:17, 21:27)]
TARGET PEATAUTO
modelo_peatauto_ls = lm(PeatAuto_std ~., df_peatauto)
summary(modelo_peatauto_ls)##
## Call:
## lm(formula = PeatAuto_std ~ ., data = df_peatauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.947 -9.491 -0.023 13.550 31.458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.2207 5.3963 12.457 3.19e-08 ***
## Poblacion 9.6839 8.8687 1.092 0.2963
## DenPob 1.2058 9.7536 0.124 0.9037
## ArCiclista -11.9114 7.2237 -1.649 0.1251
## ArBajaVel -11.3210 8.2395 -1.374 0.1946
## PMPeatones 436.1806 220.4523 1.979 0.0713 .
## PMCiclistas 102.0832 52.6049 1.941 0.0762 .
## PMTPublico 377.1615 184.6464 2.043 0.0637 .
## PMVMotor 514.9586 250.2670 2.058 0.0620 .
## Precipitacion -0.2775 8.1269 -0.034 0.9733
## Temp -8.7025 8.2750 -1.052 0.3137
## PBI -1.3527 6.7087 -0.202 0.8436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.44 on 12 degrees of freedom
## Multiple R-squared: 0.7489, Adjusted R-squared: 0.5187
## F-statistic: 3.254 on 11 and 12 DF, p-value: 0.02693Observamos un R2adj del orden del 50%, y un pvalue global, no tan bueno como nos gustaria, pero significativo. Tambien notamos que varios estimadores no resultan significativamente distintos de 0, esto puede deberse a problemas de alta colinealidad entre covariables. Vamos a realizar un analisis de la colinealidad entre las covariables.
Vamos a realizar un analisis VIF
nota: Aplicamos una transformacion logaritmica para una visualizacion mas amigable, los valores por sobre la linea roja, indican correlacion elevada.
barplot(log(vif(modelo_peatauto_ls))[vif(modelo_peatauto_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
(Cabe destacar en el gráfico se muestra el logaritmo de los valores del VIF, para facilitar la visualización)
Se observan problemas de colinealidad entre las participaciones modales de las formas de transporte: PMPeatones, PMCiclistas, PMTPublic y PMVmotor. Esto tiene logica desde la construccion de las variables modales.
Vamos a trabajar para intentar solucionar los problemas de multicolinealidad hallados.
Un primer approach podemos hacerlo a partir de metodos de seleccion de covariables como el criterio stepwise.
METODO STEPWISE TARGET PEATAUTO
modelo_peatauto_wise=step(modelo_peatauto_ls, direction = c("both"), trace = F)
summary(modelo_peatauto_wise)##
## Call:
## lm(formula = PeatAuto_std ~ Poblacion + ArCiclista + ArBajaVel +
## PMPeatones + PMCiclistas + PMTPublico + PMVMotor + Temp,
## data = df_peatauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.897 -9.821 0.269 14.257 33.062
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.221 4.836 13.899 5.66e-10 ***
## Poblacion 9.330 6.850 1.362 0.1933
## ArCiclista -12.101 6.240 -1.939 0.0715 .
## ArBajaVel -10.705 6.491 -1.649 0.1199
## PMPeatones 424.800 183.052 2.321 0.0348 *
## PMCiclistas 99.486 43.991 2.262 0.0390 *
## PMTPublico 368.083 153.363 2.400 0.0298 *
## PMVMotor 501.859 208.362 2.409 0.0293 *
## Temp -7.922 6.092 -1.300 0.2131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.69 on 15 degrees of freedom
## Multiple R-squared: 0.7479, Adjusted R-squared: 0.6135
## F-statistic: 5.563 on 8 and 15 DF, p-value: 0.002166Pasamos de tener 11 a 8 covariables en nuestro modelo, pero observamos que el metodo opto por quedarse con las covariables de participacion modal. Tambien, observamos una mejoria en el R2adj, asi como tambien en el RSE. Tambien logramos aumentar la cantidad de estimadores estadisticamente distintos de cero.
barplot(log(vif(modelo_peatauto_wise))[vif(modelo_peatauto_wise) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
Seguimos manteniendo el problema de multicolinealidad, vamos a probar otro enfoque.
Vamos a aplicar un metodo de regularizacion, para quitar covariables, a ver si logramos mejorar nuestro modelo.
METODO DE REGULARIZACION LASSO TARGET PEATAUTO
xx = model.matrix(PeatAuto_std ~., data = df_peatauto)[, -1]
modelo_peatauto_lassocv <- cv.glmnet(xx, y = df_peatauto$PeatAuto_std,
alpha = 1,
nfolds = 5)modelo_peatauto_lassocv_min <- glmnet(xx, y = df_peatauto$PeatAuto_std,
alpha = 1,
lambda = modelo_peatauto_lassocv$lambda.min)
modelo_peatauto_lassocv_1sd <- glmnet(xx, y = df_peatauto$PeatAuto_std,
alpha = 1,
lambda = modelo_peatauto_lassocv$lambda.1se)
cbind(coef(modelo_peatauto_lassocv_min),
coef(modelo_peatauto_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) 67.22071 67.220709
## Poblacion . .
## DenPob . .
## ArCiclista -11.73146 -5.340420
## ArBajaVel . .
## PMPeatones -13.59697 -7.205842
## PMCiclistas . .
## PMTPublico . .
## PMVMotor . .
## Precipitacion . .
## Temp . .
## PBI . .Ambos criterios de penalidad deciden conservar las mismas covariables para la regresion
modelo_peatauto_ls_2=lm(PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto)
summary(modelo_peatauto_ls_2)##
## Call:
## lm(formula = PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.304 -22.595 2.063 14.025 43.113
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.221 5.192 12.948 1.77e-11 ***
## ArCiclista -16.609 5.781 -2.873 0.00911 **
## PMPeatones -18.475 5.781 -3.196 0.00435 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.43 on 21 degrees of freedom
## Multiple R-squared: 0.5933, Adjusted R-squared: 0.5545
## F-statistic: 15.32 on 2 and 21 DF, p-value: 7.899e-05Observamos pvalue global e individuales mucho mas significativos. Nuestro R2adj se ubica en un 55% para el target PeatAuto_std.
par(mfrow=c(2,2))
plot(modelo_peatauto_ls_2)
Revisamos si los residuos poseen algun tipo de estructura o si se encuentra algun outlier, pero todo parece estar en orden.
control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")
modelo_MM_peatauto= lmrobdetMM(PeatAuto_std ~ ArCiclista +
PMPeatones, data = df_peatauto, control = control)
summary(modelo_MM_peatauto)##
## Call:
## lmrobdetMM(formula = PeatAuto_std ~ ArCiclista + PMPeatones, data = df_peatauto,
## control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -37.507 -20.799 2.045 15.554 43.535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.761 6.060 11.018 3.46e-10 ***
## ArCiclista -16.000 6.561 -2.438 0.0237 *
## PMPeatones -17.263 6.565 -2.630 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 30.81
## Multiple R-squared: 0.6318, Adjusted R-squared: 0.5967
## Convergence in 11 IRWLS iterationsplot(modelo_MM_peatauto$residuals, modelo_MM_peatauto$rweights,
main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")
No se observan outliers, pero conseguimos un mejor R2adj que con el metodo de LS.
TARGET CICAUTO
modelo_cicauto_ls = lm(CicAuto_std ~., data = df_cicauto)
summary(modelo_cicauto_ls)##
## Call:
## lm(formula = CicAuto_std ~ ., data = df_cicauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.807 -9.711 4.135 9.997 24.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.2005 4.4369 6.581 2.61e-05 ***
## Poblacion 3.7392 7.2920 0.513 0.617
## DenPob 0.1312 8.0196 0.016 0.987
## ArCiclista -3.2262 5.9395 -0.543 0.597
## ArBajaVel 2.1050 6.7747 0.311 0.761
## PMPeatones 233.0700 181.2602 1.286 0.223
## PMCiclistas 63.8725 43.2528 1.477 0.166
## PMTPublico 208.2031 151.8199 1.371 0.195
## PMVMotor 278.5934 205.7745 1.354 0.201
## Precipitacion -2.0798 6.6821 -0.311 0.761
## Temp -9.3761 6.8038 -1.378 0.193
## PBI -0.0853 5.5161 -0.015 0.988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.74 on 12 degrees of freedom
## Multiple R-squared: 0.6285, Adjusted R-squared: 0.288
## F-statistic: 1.846 on 11 and 12 DF, p-value: 0.1535barplot(log(vif(modelo_cicauto_ls))[vif(modelo_cicauto_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
METODO STEPWISE TARGET CICAUTO
modelo_cicauto_wise=step(modelo_cicauto_ls, direction = c("both"), trace = F)
summary(modelo_cicauto_wise)##
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + PMCiclistas + PMTPublico +
## PMVMotor + Temp, data = df_cicauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.058 -11.882 2.272 7.486 36.046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.201 3.832 7.621 4.86e-07 ***
## PMPeatones 244.590 125.882 1.943 0.0678 .
## PMCiclistas 65.834 30.005 2.194 0.0416 *
## PMTPublico 220.626 106.099 2.079 0.0521 .
## PMVMotor 291.050 144.061 2.020 0.0585 .
## Temp -6.907 4.424 -1.561 0.1359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.77 on 18 degrees of freedom
## Multiple R-squared: 0.5845, Adjusted R-squared: 0.469
## F-statistic: 5.064 on 5 and 18 DF, p-value: 0.004519METODO DE REGULARIZACION LASSO TARGET CICAUTO
xx = model.matrix(CicAuto_std ~., data = df_cicauto)[, -1]
modelo_cicauto_lassocv <- cv.glmnet(xx, y = df_cicauto$CicAuto_std,
alpha = 1,
nfolds = 5)modelo_cicauto_lassocv_min <- glmnet(xx, y = df_cicauto$CicAuto_std,
alpha = 1,
lambda = modelo_cicauto_lassocv$lambda.min)
modelo_cicauto_lassocv_1sd <- glmnet(xx, y = df_cicauto$CicAuto_std,
alpha = 1,
lambda = modelo_cicauto_lassocv$lambda.1se)
cbind(coef(modelo_cicauto_lassocv_min),
coef(modelo_cicauto_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) 29.2005120 29.20051
## Poblacion . 0.00000
## DenPob . .
## ArCiclista . .
## ArBajaVel . .
## PMPeatones -4.6539351 .
## PMCiclistas . .
## PMTPublico 0.1524345 .
## PMVMotor . .
## Precipitacion . .
## Temp -5.7587640 .
## PBI . .Observamos que el criterio 1sd de lambda, directamente nos deja sin covariables y que el criterio min, nos sugiere que dejemos tanto PMPeatones como PMTPublico, dado que ambas tienen colinealidad como marcamos previamente, optamos por dejar solo PMPeatones
modelo_cicauto_ls_2=lm(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto)
summary(modelo_cicauto_ls_2)##
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.042 -4.402 -0.553 5.602 63.282
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.201 4.210 6.935 7.5e-07 ***
## PMPeatones -9.680 4.529 -2.137 0.0445 *
## Temp -10.776 4.529 -2.379 0.0269 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.63 on 21 degrees of freedom
## Multiple R-squared: 0.4146, Adjusted R-squared: 0.3589
## F-statistic: 7.437 on 2 and 21 DF, p-value: 0.003616No nos queda un modelo muy bueno, por ahora nos quedariamos con modelo_cicauto_wise
par(mfrow=c(2,2))
plot(modelo_cicauto_wise)
Se observan registros a tener en cuenta con alto leverage y residuo cercano a 3
modelo_MM_cicauto = lmrobdetMM(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto, control = control)
summary(modelo_MM_cicauto)##
## Call:
## lmrobdetMM(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto,
## control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -35.7783 -7.2515 -0.4663 4.4901 61.3931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.090 2.961 10.500 8.2e-10 ***
## PMPeatones -12.388 3.498 -3.542 0.00193 **
## Temp -11.656 3.556 -3.278 0.00359 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 12.02
## Multiple R-squared: 0.4775, Adjusted R-squared: 0.4277
## Convergence in 68 IRWLS iterationsplot(modelo_MM_cicauto$residuals, modelo_MM_cicauto$rweights,
main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")
Encontramos 1 outlier, al que la regresion robusta le asigna un rweight de 0, y que tiene un residuo muy elevado.
sort(modelo_MM_cicauto$rweights)[1]## 20
## 0Dropeamos esta ciudad y creamos un modelo OLS:
modelo_cicauto_ls_3 = lm(CicAuto_std ~ PMPeatones + Temp, data = df_cicauto[-c(20),])
summary(modelo_cicauto_ls_3)##
## Call:
## lm(formula = CicAuto_std ~ PMPeatones + Temp, data = df_cicauto[-c(20),
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.246 -2.486 1.355 6.310 26.473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.437 3.210 8.237 7.4e-08 ***
## PMPeatones -10.205 3.381 -3.018 0.00680 **
## Temp -9.870 3.386 -2.915 0.00856 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.39 on 20 degrees of freedom
## Multiple R-squared: 0.5622, Adjusted R-squared: 0.5184
## F-statistic: 12.84 on 2 and 20 DF, p-value: 0.0002587TARGET V2RMSM
modelo_v2rmsm_ls = lm(V2RMSM_std ~., data = df_v2rmsm)
summary(modelo_v2rmsm_ls)##
## Call:
## lm(formula = V2RMSM_std ~ ., data = df_v2rmsm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3564 -5.9630 0.5517 4.1226 16.6726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.78329 7.85880 2.899 0.0145 *
## Poblacion 5.79059 3.62827 1.596 0.1388
## DenPob -0.61035 3.80811 -0.160 0.8756
## ArCiclista -0.71819 3.12160 -0.230 0.8223
## ArBajaVel -7.88648 3.45857 -2.280 0.0435 *
## PMPeatones 145.09200 99.05874 1.465 0.1710
## PMCiclistas 33.30006 23.52679 1.415 0.1846
## PMTPublico 122.81272 83.63822 1.468 0.1700
## PMVMotor 167.35331 113.57767 1.473 0.1687
## Precipitacion -3.69178 3.17114 -1.164 0.2690
## Temp 4.99689 3.37427 1.481 0.1667
## PBI 0.18098 2.62207 0.069 0.9462
## PeatAuto_std -0.09294 0.11264 -0.825 0.4268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.32 on 11 degrees of freedom
## Multiple R-squared: 0.7151, Adjusted R-squared: 0.4043
## F-statistic: 2.301 on 12 and 11 DF, p-value: 0.08913barplot(log(vif(modelo_v2rmsm_ls))[vif(modelo_v2rmsm_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
METODO STEPWISE TARGET V2RMSM
modelo_v2rmsm_wise=step(modelo_v2rmsm_ls, direction = c("both"), trace = F)
summary(modelo_v2rmsm_wise)##
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMPeatones +
## PMCiclistas + PMTPublico + PMVMotor + Precipitacion + Temp,
## data = df_v2rmsm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.6412 -6.6659 -0.4107 4.1672 16.8782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.536 1.863 8.876 2.34e-07 ***
## Poblacion 4.818 2.890 1.667 0.1162
## ArBajaVel -7.089 2.554 -2.775 0.0142 *
## PMPeatones 105.585 71.474 1.477 0.1603
## PMCiclistas 24.200 17.200 1.407 0.1798
## PMTPublico 88.430 60.056 1.472 0.1616
## PMVMotor 120.774 81.395 1.484 0.1586
## Precipitacion -3.854 2.721 -1.416 0.1771
## Temp 5.468 2.446 2.235 0.0410 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.127 on 15 degrees of freedom
## Multiple R-squared: 0.6959, Adjusted R-squared: 0.5337
## F-statistic: 4.29 on 8 and 15 DF, p-value: 0.007388METODO DE REGULARIZACION LASSO TARGET V2RMSM
xx = model.matrix(V2RMSM_std ~., data = df_v2rmsm)[, -1]
modelo_v2rmsm_lassocv <- cv.glmnet(xx, y = df_v2rmsm$V2RMSM_std,
alpha = 1,
nfolds = 5)modelo_v2rmsm_lassocv_min <- glmnet(xx, y = df_v2rmsm$V2RMSM_std,
alpha = 1,
lambda = modelo_v2rmsm_lassocv$lambda.min)
modelo_v2rmsm_lassocv_1sd <- glmnet(xx, y = df_v2rmsm$V2RMSM_std,
alpha = 1,
lambda = modelo_v2rmsm_lassocv$lambda.1se)
cbind(coef(modelo_v2rmsm_lassocv_min),
coef(modelo_v2rmsm_lassocv_1sd))## 13 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) 16.5359592 16.5359592
## Poblacion 3.1086528 0.0722265
## DenPob . .
## ArCiclista . .
## ArBajaVel -3.5556674 .
## PMPeatones . .
## PMCiclistas -0.3550315 .
## PMTPublico . .
## PMVMotor . .
## Precipitacion -1.9226689 .
## Temp 5.0761962 3.0827599
## PBI . .
## PeatAuto_std . .Vamos a adoptar el criterio min en este caso
modelo_v2rmsm_ls_2=lm(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas+
Precipitacion + Temp, data = df_v2rmsm)
summary(modelo_v2rmsm_ls_2)##
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas +
## Precipitacion + Temp, data = df_v2rmsm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5083 -6.7125 0.1061 3.5711 20.1013
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.536 1.825 9.061 3.98e-08 ***
## Poblacion 4.061 2.158 1.882 0.0761 .
## ArBajaVel -5.375 2.041 -2.634 0.0169 *
## PMCiclistas -1.039 2.034 -0.511 0.6156
## Precipitacion -3.227 2.439 -1.323 0.2023
## Temp 5.461 2.248 2.429 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.94 on 18 degrees of freedom
## Multiple R-squared: 0.6498, Adjusted R-squared: 0.5525
## F-statistic: 6.68 on 5 and 18 DF, p-value: 0.001108par(mfrow=c(2,2))
plot(modelo_v2rmsm_ls_2)
Hay algunas ciudades con alta distancia de cook y residuo superior a modulo de 2.
modelo_MM_v2rmsm = lmrobdetMM(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas +
Precipitacion + Temp, data = df_v2rmsm, control = control)
summary(modelo_MM_v2rmsm)##
## Call:
## lmrobdetMM(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas + Precipitacion +
## Temp, data = df_v2rmsm, control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -7.7589 -3.4036 -0.4829 4.6299 26.8394
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.5173 1.6220 8.950 4.78e-08 ***
## Poblacion 3.3130 2.0732 1.598 0.1275
## ArBajaVel -3.3675 1.7910 -1.880 0.0764 .
## PMCiclistas -0.8955 1.6934 -0.529 0.6034
## Precipitacion -2.5392 2.2508 -1.128 0.2741
## Temp 3.6217 1.9652 1.843 0.0819 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 8.773
## Multiple R-squared: 0.5196, Adjusted R-squared: 0.3861
## Convergence in 37 IRWLS iterationsplot(modelo_MM_v2rmsm$residuals, modelo_MM_v2rmsm$rweights,
main = "ESTUDIO DE OUTLIERS", xlab = "RESIDUOS", ylab = "RWEIGHTS", col = "blue")
sort(modelo_MM_v2rmsm$rweights)[1]## 6
## 0.001668902SegVial[c(6) ,"Ciudad"]## # A tibble: 1 × 1
## Ciudad
## <chr>
## 1 MarseilleEncontramos que la ciudad de Marseille, tiene un residuo muy alto y peso bajo segun nuestro modelo robusto. Vamos a dropearla y crear otro modelo LS, sin esta observacion.
modelo_v2rmsm_ls_3 = lm(V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas +
Precipitacion + Temp, data = df_v2rmsm[-c(6),])
summary(modelo_v2rmsm_ls_3)##
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas +
## Precipitacion + Temp, data = df_v2rmsm[-c(6), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.5434 -5.4103 0.0495 3.8454 13.9014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.4401 1.5337 10.067 1.4e-08 ***
## Poblacion 5.2117 1.8033 2.890 0.0102 *
## ArBajaVel -4.2813 1.7061 -2.509 0.0225 *
## PMCiclistas -0.4932 1.6738 -0.295 0.7718
## Precipitacion -1.0105 2.1171 -0.477 0.6392
## Temp 4.9617 1.8466 2.687 0.0156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.317 on 17 degrees of freedom
## Multiple R-squared: 0.6734, Adjusted R-squared: 0.5774
## F-statistic: 7.011 on 5 and 17 DF, p-value: 0.001008Vamos a ver si podemos encontrar alguna transformacion de las covariables que mejore aun mas nuestro modelo.
#
df_v2rmsm$PMCiclistas = df_v2rmsm$PMCiclistas + 1
df_v2rmsm$Precipitacion = df_v2rmsm$Precipitacion + 2.1
df_v2rmsm$Temp = df_v2rmsm$Temp + 2
#
# boxTidwell(V2RMSM_std ~ Precipitacion + PMCiclistas + Temp, data = df_v2rmsm[-c(6), ])df_v2rmsm$PMCiclistas2 = 1/sqrt(df_v2rmsm$PMCiclistas)
df_v2rmsm$Temp = sqrt(df_v2rmsm$Temp)
df_v2rmsm$Precipitacion2 = df_v2rmsm$Precipitacion^2
modelo_v2rmsm_ls_3 = lm(V2RMSM_std ~ Poblacion + ArBajaVel +
PMCiclistas2 + Precipitacion2 + Temp , data = df_v2rmsm[-c(6),])
summary(modelo_v2rmsm_ls_3)##
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + PMCiclistas2 +
## Precipitacion2 + Temp, data = df_v2rmsm[-c(6), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6600 -5.3257 -0.1406 4.0773 13.1634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1745 7.9853 0.272 0.78866
## Poblacion 5.7118 1.7797 3.209 0.00514 **
## ArBajaVel -5.3317 1.8920 -2.818 0.01185 *
## PMCiclistas2 -2.3264 4.8725 -0.477 0.63912
## Precipitacion2 -0.3877 0.3343 -1.160 0.26225
## Temp 13.3211 4.4019 3.026 0.00762 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.943 on 17 degrees of freedom
## Multiple R-squared: 0.7059, Adjusted R-squared: 0.6195
## F-statistic: 8.163 on 5 and 17 DF, p-value: 0.000439modelo_v2rmsm_wise2 = step(modelo_v2rmsm_ls_3, direction = c("both"), trace = F)
summary(modelo_v2rmsm_wise2)##
## Call:
## lm(formula = V2RMSM_std ~ Poblacion + ArBajaVel + Precipitacion2 +
## Temp, data = df_v2rmsm[-c(6), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.7252 -5.1752 0.3439 3.2655 13.7347
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1084 6.5656 0.017 0.98701
## Poblacion 5.3045 1.5281 3.471 0.00272 **
## ArBajaVel -4.8303 1.5396 -3.137 0.00569 **
## Precipitacion2 -0.4162 0.3218 -1.293 0.21220
## Temp 12.8700 4.2061 3.060 0.00675 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.792 on 18 degrees of freedom
## Multiple R-squared: 0.702, Adjusted R-squared: 0.6358
## F-statistic: 10.6 on 4 and 18 DF, p-value: 0.0001356TARGET V2RMAUTO
modelo_v2rmauto_ls = lm(V2RMAuto_std ~., data = df_v2rmauto)
summary(modelo_v2rmauto_ls)##
## Call:
## lm(formula = V2RMAuto_std ~ ., data = df_v2rmauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.083 -14.225 -0.007 12.132 28.776
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.263 4.931 12.627 2.74e-08 ***
## Poblacion 13.992 8.104 1.727 0.1099
## DenPob -5.643 8.913 -0.633 0.5385
## ArCiclista -13.765 6.601 -2.085 0.0591 .
## ArBajaVel -12.994 7.529 -1.726 0.1100
## PMPeatones 514.017 201.448 2.552 0.0254 *
## PMCiclistas 121.364 48.070 2.525 0.0267 *
## PMTPublico 418.179 168.728 2.478 0.0290 *
## PMVMotor 578.256 228.692 2.529 0.0265 *
## Precipitacion -7.792 7.426 -1.049 0.3148
## Temp 20.817 7.562 2.753 0.0175 *
## PBI -2.411 6.130 -0.393 0.7010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.16 on 12 degrees of freedom
## Multiple R-squared: 0.8031, Adjusted R-squared: 0.6227
## F-statistic: 4.451 on 11 and 12 DF, p-value: 0.008055modelo_v2rmauto_wise = step(modelo_v2rmauto_ls, direction = c("both"),
trace = F)
summary(modelo_v2rmauto_wise)##
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + ArCiclista + ArBajaVel +
## PMPeatones + PMCiclistas + PMTPublico + PMVMotor + Precipitacion +
## Temp, data = df_v2rmauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.676 -15.748 -2.734 10.980 32.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.263 4.693 13.266 2.55e-09 ***
## Poblacion 12.570 7.448 1.688 0.11361
## ArCiclista -14.850 6.127 -2.424 0.02950 *
## ArBajaVel -14.393 6.456 -2.229 0.04268 *
## PMPeatones 469.573 180.098 2.607 0.02068 *
## PMCiclistas 111.695 43.346 2.577 0.02195 *
## PMTPublico 381.361 151.355 2.520 0.02452 *
## PMVMotor 530.756 205.143 2.587 0.02150 *
## Precipitacion -8.022 6.936 -1.156 0.26684
## Temp 19.378 6.163 3.144 0.00717 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.99 on 14 degrees of freedom
## Multiple R-squared: 0.7919, Adjusted R-squared: 0.6582
## F-statistic: 5.921 on 9 and 14 DF, p-value: 0.001701barplot(log(vif(modelo_v2rmauto_wise))[vif(modelo_v2rmauto_wise) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
METODO DE REGULARIZACION LASSO TARGET V2RMAUTO
xx = model.matrix(V2RMAuto_std ~., data = df_v2rmauto)[, -1]
modelo_v2rmauto_lassocv <- cv.glmnet(xx, y = df_v2rmauto$V2RMAuto_std,
alpha = 1,
nfolds = 5)modelo_v2rmauto_lassocv_min <- glmnet(xx, y = df_v2rmauto$V2RMAuto_std,
alpha = 1,
lambda = modelo_v2rmauto_lassocv$lambda.min)
modelo_v2rmauto_lassocv_1sd <- glmnet(xx, y = df_v2rmauto$V2RMAuto_std,
alpha = 1,
lambda = modelo_v2rmauto_lassocv$lambda.1se)
cbind(coef(modelo_v2rmauto_lassocv_min),
coef(modelo_v2rmauto_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) 62.2628935 62.26289
## Poblacion 0.6046567 .
## DenPob . .
## ArCiclista -7.9051384 .
## ArBajaVel . .
## PMPeatones . .
## PMCiclistas . .
## PMTPublico . .
## PMVMotor . .
## Precipitacion -2.7152213 .
## Temp 19.6160028 10.04084
## PBI . .Vamos a adoptar el criterio min en este caso
modelo_v2rmauto_ls_2=lm(V2RMAuto_std ~ Poblacion + ArCiclista +
Precipitacion + Temp, data = df_v2rmauto)
summary(modelo_v2rmauto_ls_2)##
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + ArCiclista + Precipitacion +
## Temp, data = df_v2rmauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.942 -9.842 -1.198 11.807 54.438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.263 5.354 11.630 4.39e-10 ***
## Poblacion 2.223 6.448 0.345 0.73411
## ArCiclista -14.727 5.897 -2.498 0.02185 *
## Precipitacion -7.668 7.201 -1.065 0.30032
## Temp 22.959 6.397 3.589 0.00196 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.23 on 19 degrees of freedom
## Multiple R-squared: 0.6326, Adjusted R-squared: 0.5552
## F-statistic: 8.178 on 4 and 19 DF, p-value: 0.0005183modelo_v2rmauto_wise_2=step(modelo_v2rmauto_ls_2, direction = c("both"),
trace = F)
summary(modelo_v2rmauto_wise_2)##
## Call:
## lm(formula = V2RMAuto_std ~ ArCiclista + Precipitacion + Temp,
## data = df_v2rmauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.636 -11.680 -0.996 11.353 51.952
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.263 5.234 11.895 1.59e-10 ***
## ArCiclista -15.379 5.461 -2.816 0.01067 *
## Precipitacion -8.731 6.362 -1.372 0.18515
## Temp 23.043 6.250 3.687 0.00146 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.64 on 20 degrees of freedom
## Multiple R-squared: 0.6303, Adjusted R-squared: 0.5748
## F-statistic: 11.37 on 3 and 20 DF, p-value: 0.0001438modelo_MM_v2rmauto= lmrobdetMM(V2RMAuto_std ~ Poblacion + ArCiclista + Precipitacion +
Temp, data = df_v2rmauto, control = control)
plot(modelo_MM_v2rmauto$residuals, modelo_MM_v2rmauto$rweights)
Observamos dos ciudades que son outliers
sort(modelo_MM_v2rmauto$rweights)[1:2]## 6 9
## 0 0modelo_v2rmauto_wise_3=step(lm(V2RMAuto_std ~ ., data = df_v2rmauto[-c(6, 9),]),
direction = c("both"),
trace = F)
summary(modelo_v2rmauto_wise_3)##
## Call:
## lm(formula = V2RMAuto_std ~ Poblacion + PMVMotor + Temp + PBI,
## data = df_v2rmauto[-c(6, 9), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.094 -6.085 -1.078 7.216 32.995
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 54.166 2.853 18.985 7.00e-13 ***
## Poblacion 18.051 3.257 5.543 3.58e-05 ***
## PMVMotor 10.612 3.108 3.415 0.00330 **
## Temp 11.825 3.242 3.647 0.00199 **
## PBI 3.958 2.918 1.356 0.19275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.21 on 17 degrees of freedom
## Multiple R-squared: 0.8061, Adjusted R-squared: 0.7605
## F-statistic: 17.67 on 4 and 17 DF, p-value: 6.9e-06par(mfrow=c(2,2))
plot(modelo_v2rmauto_wise_3)
No se observa estructura, ni residuos mayores a modulos de 3.
TARGET AUTOSM
modelo_autosm_ls=lm(AutoSM_std ~., data = df_autosm)
summary(modelo_autosm_ls)##
## Call:
## lm(formula = AutoSM_std ~ ., data = df_autosm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.7315 -4.7622 0.5719 5.7051 16.7060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.6636 2.1624 8.169 3.03e-06 ***
## Poblacion -0.5314 3.5538 -0.150 0.8836
## DenPob -3.7350 3.9084 -0.956 0.3581
## ArCiclista -3.0996 2.8946 -1.071 0.3053
## ArBajaVel -5.9694 3.3017 -1.808 0.0957 .
## PMPeatones 49.8306 88.3383 0.564 0.5831
## PMCiclistas 12.6611 21.0795 0.601 0.5593
## PMTPublico 47.7785 73.9904 0.646 0.5306
## PMVMotor 61.0973 100.2855 0.609 0.5537
## Precipitacion -2.4117 3.2566 -0.741 0.4732
## Temp 3.0742 3.3159 0.927 0.3721
## PBI 2.1197 2.6883 0.788 0.4457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.59 on 12 degrees of freedom
## Multiple R-squared: 0.663, Adjusted R-squared: 0.3541
## F-statistic: 2.146 on 11 and 12 DF, p-value: 0.1026barplot(log(vif(modelo_autosm_ls))[vif(modelo_autosm_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
METODO DE REGULARIZACION LASSO TARGET AUTOSM
xx = model.matrix(AutoSM_std ~., data = df_autosm)[, -1]
modelo_autosm_lassocv <- cv.glmnet(xx, y = df_autosm$AutoSM_std,
alpha = 1,
nfolds = 5,
maxit=1000000)modelo_autosm_lassocv_min <- glmnet(xx, y = df_autosm$AutoSM_std,
alpha = 1,
lambda = modelo_autosm_lassocv$lambda.min)
modelo_autosm_lassocv_1sd <- glmnet(xx, y = df_autosm$AutoSM_std,
alpha = 1,
lambda = modelo_autosm_lassocv$lambda.1se)
cbind(coef(modelo_autosm_lassocv_min),
coef(modelo_autosm_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) 17.663610 17.6636102
## Poblacion . .
## DenPob . .
## ArCiclista -1.763633 -0.5199541
## ArBajaVel -4.210453 -2.9179391
## PMPeatones -2.254634 -1.2171255
## PMCiclistas . .
## PMTPublico . .
## PMVMotor . .
## Precipitacion . .
## Temp . .
## PBI . .modelo_autosm_ls2 = lm(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones,
data = df_autosm)
summary(modelo_autosm_ls2)##
## Call:
## lm(formula = AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones,
## data = df_autosm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.927 -6.676 -1.188 6.273 17.445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.664 1.919 9.203 1.25e-08 ***
## ArCiclista -3.427 2.160 -1.587 0.1282
## ArBajaVel -5.940 2.127 -2.793 0.0112 *
## PMPeatones -3.643 2.234 -1.631 0.1186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.402 on 20 degrees of freedom
## Multiple R-squared: 0.5575, Adjusted R-squared: 0.4912
## F-statistic: 8.401 on 3 and 20 DF, p-value: 0.0008222par(mfrow=c(2,2))
plot(modelo_autosm_ls2)
A partir de los graficos de los residuos, decidimos intentar alguna transformacion a ver si mejora nuestro modelo.
df_autosm$AutoSM_std = df_autosm$AutoSM_std + 0.01
boxcox(modelo_autosm_ls2)
df_autosm$AutoSM_std = sqrt(df_autosm$AutoSM_std)
modelo_autosm_ls3 = lm(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones,
data = df_autosm)
summary(modelo_autosm_ls3)##
## Call:
## lm(formula = AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones,
## data = df_autosm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.86367 -0.74229 0.06008 1.03209 1.52277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8928 0.2300 16.922 2.55e-13 ***
## ArCiclista -0.4568 0.2589 -1.765 0.09290 .
## ArBajaVel -0.7946 0.2549 -3.117 0.00543 **
## PMPeatones -0.3729 0.2677 -1.393 0.17896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.127 on 20 degrees of freedom
## Multiple R-squared: 0.5799, Adjusted R-squared: 0.5169
## F-statistic: 9.203 on 3 and 20 DF, p-value: 0.0004977par(mfrow=c(2,2))
plot(modelo_autosm_ls3)
modelo_MM_autosm = lmrobdetMM(AutoSM_std ~ ArCiclista + ArBajaVel + PMPeatones,
data = df_autosm, control = control)
plot(modelo_MM_autosm$residuals, modelo_MM_autosm$rweights)
No vemos outliers a considerar
TARGET AUTOAUTO
modelo_autoauto_ls = lm(AutoAuto_std ~., data = df_autoauto)
summary(modelo_autoauto_ls)##
## Call:
## lm(formula = AutoAuto_std ~ ., data = df_autoauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99348 -0.37553 0.00617 0.24130 2.51788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.137e-15 1.958e-01 0.000 1.000
## Poblacion -1.467e-01 3.218e-01 -0.456 0.657
## DenPob 1.595e-01 3.539e-01 0.451 0.660
## ArCiclista -1.533e-01 2.621e-01 -0.585 0.570
## ArBajaVel -2.961e-01 2.990e-01 -0.990 0.342
## PMPeatones 1.215e+00 7.999e+00 0.152 0.882
## PMCiclistas 2.535e-01 1.909e+00 0.133 0.897
## PMTPublico 1.536e+00 6.700e+00 0.229 0.823
## PMVMotor 1.965e+00 9.081e+00 0.216 0.832
## Precipitacion -3.464e-01 2.949e-01 -1.175 0.263
## Temp -2.556e-01 3.003e-01 -0.851 0.411
## PBI -7.630e-02 2.434e-01 -0.313 0.759
##
## Residual standard error: 0.9593 on 12 degrees of freedom
## Multiple R-squared: 0.5199, Adjusted R-squared: 0.07983
## F-statistic: 1.181 on 11 and 12 DF, p-value: 0.3878Realizamos un estudio de la colinealidad de las covariables del modelo, ya que tan pocos betas significativos nos da la idea de un problema por ese lado.
barplot(log(vif(modelo_autoauto_ls))[vif(modelo_autoauto_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
APLICAMOS METODO LASSO DE REGULARIZACION TARGET AUTOAUTO
xx = model.matrix(AutoAuto_std ~., data = df_autoauto)[, -1]
modelo_autoauto_lassocv <- cv.glmnet(xx, y = df_autoauto$AutoAuto_std,
alpha = 1,
nfolds = 5)modelo_autoauto_lassocv_min <- glmnet(xx, y = df_autoauto$AutoAuto_std,
alpha = 1,
lambda = modelo_autoauto_lassocv$lambda.min)
modelo_autoauto_lassocv_1sd <- glmnet(xx, y = df_autoauto$AutoAuto_std,
alpha = 1,
lambda = modelo_autoauto_lassocv$lambda.1se)
cbind(coef(modelo_autoauto_lassocv_min),
coef(modelo_autoauto_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) -1.698226e-16 -1.457168e-16
## Poblacion . 0.000000e+00
## DenPob . .
## ArCiclista . .
## ArBajaVel . .
## PMPeatones -2.895011e-01 .
## PMCiclistas . .
## PMTPublico . .
## PMVMotor . .
## Precipitacion . .
## Temp . .
## PBI . .Si bien nuestro modelo lasso nos devuelve PMCiclistas y PMPeatones, vamos a optar por quedarnos solo con PMPeatones, ya que por construccion sabemos que ambas covariables tienen problemas de multicolinealidad elevada.
modelo_autoauto_ls_2 = lm(AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones, data = df_autoauto)
summary(modelo_autoauto_ls_2)##
## Call:
## lm(formula = AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones,
## data = df_autoauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0828 -0.3846 -0.2230 0.2706 2.8322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.919e-16 1.650e-01 0.000 1.0000
## ArCiclista -1.956e-01 1.857e-01 -1.053 0.3047
## ArBajaVel -1.712e-01 1.828e-01 -0.936 0.3602
## PMPeatones -4.647e-01 1.921e-01 -2.420 0.0252 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8084 on 20 degrees of freedom
## Multiple R-squared: 0.4318, Adjusted R-squared: 0.3465
## F-statistic: 5.066 on 3 and 20 DF, p-value: 0.009019summary(step(modelo_autoauto_ls_2, direction = c("both"), trace = F))##
## Call:
## lm(formula = AutoAuto_std ~ PMPeatones, data = df_autoauto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0520 -0.5155 -0.2052 0.1780 2.9979
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.394e-16 1.662e-01 0.000 1.00000
## PMPeatones -6.049e-01 1.698e-01 -3.563 0.00174 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8142 on 22 degrees of freedom
## Multiple R-squared: 0.3659, Adjusted R-squared: 0.3371
## F-statistic: 12.69 on 1 and 22 DF, p-value: 0.00174modelo_MM_autoauto = lmrobdetMM(AutoAuto_std ~ ArCiclista + ArBajaVel +
PMPeatones, data = df_autoauto, control = control)
plot(modelo_MM_autoauto$residuals, modelo_MM_autoauto$rweights)
Observamos 2 ciudades con valores de residuos muy elevados
sort(modelo_MM_autoauto$rweights)[1:2]## 22 13
## 0.00000000 0.01822396modelo_autoauto_ls_3=lm(AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones,
data = df_autoauto[-c(13, 22),])
summary(modelo_autoauto_ls_3)##
## Call:
## lm(formula = AutoAuto_std ~ ArCiclista + ArBajaVel + PMPeatones,
## data = df_autoauto[-c(13, 22), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81049 -0.17920 0.05436 0.21974 0.60267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.19813 0.07612 -2.603 0.01800 *
## ArCiclista -0.15063 0.08161 -1.846 0.08145 .
## ArBajaVel -0.07311 0.08097 -0.903 0.37848
## PMPeatones -0.27977 0.08756 -3.195 0.00502 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3547 on 18 degrees of freedom
## Multiple R-squared: 0.5857, Adjusted R-squared: 0.5167
## F-statistic: 8.483 on 3 and 18 DF, p-value: 0.001003par(mfrow=c(2,2))
plot(modelo_autoauto_ls_3)
Observamos una buena apariencia en los residuos
cor.test(abs(modelo_autoauto_ls_3$residuals),
modelo_autoauto_ls_3$fitted.values, method = "spearman")##
## Spearman's rank correlation rho
##
## data: abs(modelo_autoauto_ls_3$residuals) and modelo_autoauto_ls_3$fitted.values
## S = 1532, p-value = 0.5478
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.134952Comprobamos que no hay Heterocedasticidad
vif(modelo_autoauto_ls_3)## ArCiclista ArBajaVel PMPeatones
## 1.167712 1.128402 1.224928Eliminamos el problema de multicolinealidad entre las covariables
Vamos ahora a analizar todas las clases de accidentes viales, utilizando la variable accidentes_viales, que es una combinacion lineal de las anteriores
TARGET ACCIDENTES VIALES
modelo_accidentes_viales_ls = lm(accidentes_viales_std ~. ,df_accidentesviales)
summary(modelo_accidentes_viales_ls)##
## Call:
## lm(formula = accidentes_viales_std ~ ., data = df_accidentesviales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.18745 -0.44059 0.06232 0.43035 1.38822
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.477e-15 1.582e-01 0.000 1.0000
## Poblacion 2.280e-01 2.600e-01 0.877 0.3977
## DenPob -1.974e-02 2.859e-01 -0.069 0.9461
## ArCiclista -3.354e-01 2.118e-01 -1.584 0.1393
## ArBajaVel -4.172e-01 2.416e-01 -1.727 0.1097
## PMPeatones 1.228e+01 6.463e+00 1.900 0.0817 .
## PMCiclistas 2.958e+00 1.542e+00 1.918 0.0792 .
## PMTPublico 1.064e+01 5.413e+00 1.965 0.0730 .
## PMVMotor 1.445e+01 7.337e+00 1.970 0.0724 .
## Precipitacion -2.692e-01 2.383e-01 -1.130 0.2807
## Temp 1.019e-02 2.426e-01 0.042 0.9672
## PBI -4.025e-02 1.967e-01 -0.205 0.8413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.775 on 12 degrees of freedom
## Multiple R-squared: 0.6866, Adjusted R-squared: 0.3993
## F-statistic: 2.39 on 11 and 12 DF, p-value: 0.07493barplot(log(vif(modelo_accidentes_viales_ls))[vif(modelo_accidentes_viales_ls) > 5],
main = "log VIF DE LAS COVARIABLES", col = "blue")
abline(h = log(5), lwd = 2, col = "red")
APLICAMOS METODO LASSO DE REGULARIZACION TARGET ACCIDENTES_VIALES
xx = model.matrix(accidentes_viales_std ~., data = df_accidentesviales)[, -1]
modelo_accidentesviales_lassocv <- cv.glmnet(xx, y = df_accidentesviales$accidentes_viales_std,
alpha = 1,
nfolds = 5)modelo_accidentesviales_lassocv_min <- glmnet(xx, y = df_accidentesviales$accidentes_viales_std,
alpha = 1,
lambda = modelo_accidentesviales_lassocv$lambda.min)
modelo_accidentesviales_lassocv_1sd <- glmnet(xx, y = df_accidentesviales$accidentes_viales_std,
alpha = 1,
lambda = modelo_accidentesviales_lassocv$lambda.1se)
cbind(coef(modelo_accidentesviales_lassocv_min),
coef(modelo_accidentesviales_lassocv_1sd))## 12 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) -1.901963e-16 -1.656208e-16
## Poblacion 1.906163e-02 .
## DenPob . .
## ArCiclista -3.022070e-01 -1.261378e-01
## ArBajaVel -2.195146e-02 .
## PMPeatones -2.311678e-01 -5.059179e-02
## PMCiclistas . .
## PMTPublico . .
## PMVMotor . .
## Precipitacion . .
## Temp . .
## PBI . .modelo_accidentes_viales_ls2 = lm(accidentes_viales_std ~ Poblacion + ArCiclista
+ ArBajaVel + PMPeatones, data = df_accidentesviales)
summary(modelo_accidentes_viales_ls2)##
## Call:
## lm(formula = accidentes_viales_std ~ Poblacion + ArCiclista +
## ArBajaVel + PMPeatones, data = df_accidentesviales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.31256 -0.50546 -0.03908 0.40300 1.37683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.846e-17 1.505e-01 0.000 1.0000
## Poblacion 2.466e-01 1.631e-01 1.512 0.1471
## ArCiclista -3.563e-01 1.762e-01 -2.022 0.0575 .
## ArBajaVel -2.046e-01 1.715e-01 -1.193 0.2474
## PMPeatones -3.516e-01 1.756e-01 -2.002 0.0598 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7371 on 19 degrees of freedom
## Multiple R-squared: 0.5511, Adjusted R-squared: 0.4566
## F-statistic: 5.832 on 4 and 19 DF, p-value: 0.00309Observamos que los coeficientes estimados de nuestra regresion tienen un pvalue, cercano al 5%, deberemos tener en cuenta esto.
vif(modelo_accidentes_viales_ls2)## Poblacion ArCiclista ArBajaVel PMPeatones
## 1.126095 1.314217 1.244303 1.305666No se observa problemas de multicolinealidad entre las covariables.
par(mfrow=c(2,2))
plot(modelo_accidentes_viales_ls2)
Observamos que la observacion 20, tiene un muy alto leverage y un residuo de 2, es una ciudad que se debe estudiar con cuidado. Fuera de eso, no se ve estructura en los residuos y parece haber homocedasticidad y distribucion normal de los residuos. En el QQplot, volvemos a ver resaltada la observacion 22.
modelo_MM_accidentesviales = lmrobdetMM(accidentes_viales_std ~ Poblacion + ArCiclista +
ArBajaVel + PMPeatones, data = df_accidentesviales, control = control)
plot(modelo_MM_accidentesviales$residuals, modelo_MM_accidentesviales$rweights)
No se encuentran outliers. No nos encontramos que se le haya asignado un rweight de 0 a la observacion 22, con lo que optamos por no considerarla outlier y la conservaremos.
summary(modelo_MM_accidentesviales)##
## Call:
## lmrobdetMM(formula = accidentes_viales_std ~ Poblacion + ArCiclista + ArBajaVel +
## PMPeatones, data = df_accidentesviales, control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -1.334816 -0.473743 0.008481 0.402913 1.446983
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03804 0.17184 -0.221 0.827
## Poblacion 0.21387 0.19004 1.125 0.274
## ArCiclista -0.30853 0.19437 -1.587 0.129
## ArBajaVel -0.22743 0.19878 -1.144 0.267
## PMPeatones -0.32415 0.19935 -1.626 0.120
##
## Robust residual standard error: 0.8036
## Multiple R-squared: 0.5982, Adjusted R-squared: 0.5136
## Convergence in 37 IRWLS iterationsSi bien el R2adj de nuestro modelo robusto es superior al de OLS anteriormente mencionado (modelo_accidentes_viales_ls2), vamos a preferir el anterior para tareas de inferencia, debido a que sus coeficientes son mas significativos que este caso.